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NOMENCLATURE 


[Cl interface compatibility coefficient matrix 
[D] force coefficient matrix 

{F} force vector 

{Gj} term j in expansion of interface force vector 

[K] stiffness matrix 

[Ml mass matrix 

{q} modal coordinate vector 

{q} modal coordinate vector due to applied forces 

{q} modal coordinate vector due to interface forces 

t time 

ti time at time increment i 

{X} physical degrees of freedom 

At length of integration time step 

£ percent critical modal damping 

[(f)] mode shape matrix 

[wl diagonal matrix of circular frequencies 

[or] diagonal matrix of eigenvalues 

Subscripts 

a,b substructure a,b, respectively 

Superscripts 
I interface 


v 



TECHNICAL PAPER 


A TRANSIENT RESPONSE METHOD FOR LINEAR COUPLED SUBSTRUCTURES 

INTRODUCTION 


The successful design of a modern aerospace structure requires an accurate determination of the 
internal loads it will experience during its useful lifetime. Often the internal loads that impact the design 
the most are those caused by the transient response of the structure due to suddenly applied external 
forces. For example, the ignition and thrust buildup of launch vehicle engines and the subsequent lift-off 
of the vehicle results in a transient response that can cause the maximum internal loads the vehicle struc- 
ture will experience during its lifetime. Transient response analyses of structures has been the subject of 
study for many years because of its effect on the design of structures. 

The equations of motion representing a linear structural system are a set of coupled second-order 
differential equations. Even though a number of numerical methods are available for the direct solution 
of such equations, it is often not practical because of the large number of equations involved. Models of 
structures generated by finite-element computer codes often can have several thousand degrees of free- 
dom (DOF). The standard approach to overcome the problem of too many DOF is to first compute the 
vibration modes and frequencies of the structure by solving an eigenvalue problem. The structural vibra- 
tion modes are then used to transform from the original discrete physical coordinates to a set of new 
modal coordinates. This permits the number of equations to be reduced by truncating the higher fre- 
quency modal coordinates without significantly reducing the accuracy of the results. Also, the trans- 
formation results in a set of uncoupled rather than a coupled set of equations. This approach yields a set 
of equations of form and size that are easily solved on present day computers. The difficulty of perform- 
ing a vibration analysis to determine the modes and frequencies of the structure prior to solving the set of 
equations expressed in modal coordinates is still a major problem due to the large order of the eigenvalue 
problem that must be solved. 

The vibration analysis of structures has been the subject of many papers. Two approaches have 
been used to make the vibration analysis of large structural models tractable. One approach is to develop 
efficient algorithms for the solution of large eigenvalue problems directly, and the other approach is to 
make judicious approximations to reduce the size of the eigenvalue problem without significantly affect- 
ing the results. One example of the first is based on the repeated application of the Rayleigh-Ritz method 
of vibration analysis in which each solution is an improvement over the previous. This method was 
originally called an iterative Rayleigh-Ritz technique fl] and later renamed subspace iteration method 
[2,3], The other approach, called substructure synthesis, results in a reduction in size of the eigenvalue 
problem. Many vibration analysis methods (4-1 1] are based on this approach. Substructure synthesis is a 
method where the structure is considered to be made up of a number of substructures. The motion of each 
substructure is represented by a set of displacement vectors. The vectors may or may not be substructure 
modes. The number of vectors required to represent the motion of a substructure is much less than the 
number of physical DOF’s required to represent its motion. As a result of this approximation, the size of 
the eigenvalue problem is greatly reduced. 


Other transient response methods have been studied that circumvent the need to solve a large 
eigenvalue problem and still reduce the number of DOF in the equations of motion to a manageable 
number. These methods [12-14] are based on the direct integration of the equations of motion expressed 
in coordinates other than system discrete coordinates or system modal coordinates. They use a set of 
coordinates that contains both interface discrete coordinates and substructure modal coordinates. 

In this paper, a method is developed for the transient response analysis of structures that is 
different from the approaches described above. The structure is considered to be composed of substruc- 
tures. The method is based on approximating the interface forces between the substructures as power 
series in time. The equations of motion of each substructure are integrated for a time step, with all 
external forces applied including the interface forces with unknown coefficients. The unknown coeffi- 
cients in the power series are evaluated by satisfying the interface compatibility relationships at the end of 
the time step. Once the coefficients in the power series are known, the interface forces and the response 
characteristics at the end of the time step can be computed. This procedure is repeated as many times as 
required to span the time interval of interest. The Implementation of the method on a computer is quite 
simple and permits enforcing the interface compatibility relationships by a simple matrix multiplication. 
Satisfying boundary conditions between substructures in this manner allows the boundary conditions to 
be changed during the transient response by simply changing the matrix used to enforce the compatibility 
relationships. Therefore, this method is ideally suited for the transient response analysis of structures that 
have interface boundary conditions between substructures that change during the response. The transient 
response analysis of launch vehicles during the lift-off event is an example of this type of problem, since 
the boundary conditions between the vehicle/pad change as the vehicle lifts off. 


GOVERNING EQUATIONS 


For clarity, the governing equations are formulated for a structural system composed of two 
substructures. However, the extension of the equations for systems consisting of an arbitrary number of 
substructures or bodies is evident. The derivation of the equations begins by writing the equations of 
motion of the two bodies and the associated interface compatibility equations for the system shown in 
Figure I. The interface compatibility conditions for the two bodies are expressed as 


{XL} = {X’} and {Fj,} + {F' b } = {0} . (1) 

For simplicity the interface force compatibility will be rewritten as 

{F 1 ,} = {F 1 } and {F{,} - -{F 1 } . (2) 


The equations of motion for the substructures were derived by using Lagrange’s equation with 
undetermined multipliers [15], The undamped equations of motion of the substructures are written as 
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[MJ {X a } + [KJ {X ;i } = {F ; ,(t)} + 


{ 0 } 


(3a) 


[M b ] {X b } + [K b ] {X b } = {F b (t)} 



(3b) 


The interface forces are the undetermined multipliers. The equations of motion expressed by 
equations (3a) and (3b) are in discrete coordinates with the interface coordinates last. The equations can 
be solved in this form, but for large systems the computer storage requirement can be excessive. This 
problem is overcome by using the substructures vibration modes to transform the equations of motion 
from discrete coordinates to normal mode coordinates. This transformation uncouples the equations and 
results in a system of equations much simpler to solve from a computational point of view. Modal damp- 
ing of the substructures can be introduced at this point. Equations (3a) and (3b) transformed to modal 
coordinates are 


{q a } + 2Uo) a ] {q a } + K 2 ] {q a } = [D a ] {F a (t)} + [Dj] {F'(t)> , (4a) 

{qj + 2£ b [io b ] {q b } + Leo 2 ] {q b } = [D b ] {F b (t)| + [Dj,] |F T (t)} , (4b) 

where 

{X a } = [4> a l {qj and {X b } = [4> b ] {q b } . (5) 


Since the interface forces {F‘(t)} result from the motion of a structure that is governed by second order 
differential equations, they will be of trigonometric form. If equations (4a) and (4b) are solved in a step- 
wise fashion, it is appropriate to approximate the interface forces by a power series valid for a time step 
of length At. This is expressed as 


00 


{F'(t)J = 2 {Gj} (t - t,)i 


j = 0 


t; sc t *£ tj + At 


( 6 ) 


This series is expected to converge very rapidly lor time steps of the size normally used in integra- 
tion of equations of motion. Therefore, very few terms need to be retained in the series to achieve a good 
approximation ol the interlace forces. Terms up through the third order will be retained in this develop- 
ment. Truncating equation (6) and substituting in equations (4a) and (4b) yields 


{qj + 2£ a K] {q,} + K 2 ] {q a } = [D a l {F a (t)} + [D[] {{G 0 } + {G,}(t-t,) 
+ {G 2 } o-t,) 2 + {G 3 } (t-tj)- 3 } ; t, ^ t *5 t, + At , 


(7a) 


{qj + 2( b [w b j {q b } + [iog| {q b } = [D b j {F b (t)| - [D' b | {{G u } + {G,}(t-t,) 

+ {G 2 } (t-ti) 2 + {G 3 } (t-tj) 3 } ; t; *£ t tj + At . (7b) 


Since equations (7a) and (7b) are linear differential equations, superposition of their solutions are per- 
mitted. Therefore, let us define the following 


{qa} - {q a } + ft a } ; {qj = {q a } + fta} ; {q a } = {q a } + ft a } , 


{qj = {qJ + {qj ; ft b } = ft b } + ftb} ; ft h } = {q b } + ftb} , 


(8a) 

(8b) 


where {q a }, {q a }, ftbK fttl, and their derivatives are obtained by solving the following equations for tj t 
=£ tj + At 


fta} + 2C a [to a ] {qj + [u> 2 ] {q a > = [D ;i ] {F a (t)} + [D' a j {G () } , (9a) 

ft.} + 2 £ a Kl ft a } + [o> 2 ] {qj = [Di] {{G,} (t-tj) + {G 2 } (t-tj 2 + {G.,} (t-t,) 3 } , (9b) 

{qj + 2£ b [a> b ] {4,} + f a*] {q b } = [D b ] {F b (t)} - [D' b ] {G 0 } , (9c) 

ftj + 2( b [a)„] ft b } + [oo 2 b ] {q b } = [D b J {{G,} (t-t,) + {G 2 } (t-tj) 2 + {G 3 } (t-t;) 3 } , (9d) 

where the initial conditions for equations (9a) through (9d) are 
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{q a (ti)} = {q.(ti)} . {q a (ti)} = Waft)} - 

(10a) 

II 

II 

?<T 

(10b) 

{q b (ti)} = {qb(tj)} . {qb(ti)} - to)} . 

(10c) 

{q b (t,)} = {0} , to)} = i°} ’ 

(10d) 

and {G„}, the first time in the power series, is determined by substituting t, in equation (6). 

{Go} = {F T (ti)> . 

(11) 

Equations (9a) and (9c) can also be solved numerically ■ 
tj + , = tj + At. Closed-form solutions are usually possible since 
are often simple expressions. These solutions are 

or closed form to obtain the solutions at 
the applied forces over a time increment 

{qadi+i)} ; toi+i)} ; to, .)} » 

(12a) 

{q b (ti , ,)} ; {q b (h , i)} » Iqbdi , i )} 

(12b) 

Equations (9b) and (9d) can be solved in closed form, however, the solutions contain {G,}, {G 2 }, 
and {G 3 } as unknowns. The use of superposition of solutions of linear differential equations is the most 
practical approach for getting the solutions in a form suitable for evaluating the unknowns {G,}, {G 2 }, and 
{G 3 }. This is done by solving equations (9b) and (9d), with the elements in {G|}, {G 2 }, and {G 3 } being 
assigned a unit value, one at a time, and summing the solutions. This can be expressed in matrix torm as 

ffia(t i+ ,)} = [Call {G,} + [C a2 ] {G 2 } + [C a3 ] {G 3 } 

(13a) 


{qfaOi+i)} = [C ;1 |] {Gj + [C a2 l {G 2 } + [C a3 i {G 3 } , (13b) 

{g a (t i+ ,)} = [C al ] {G,} + [C a2 ] {G 2 } + [C a3 ] {G 3 } 


(13c) 


ffibdi, ■)} = [c b ,i {G,} + |C b2 l {G 2 } + [C M ] {g 3 } 


(14a) 


{«|b(ti + .)} = [C b i] {G,} + [C b2 ] {G 2 } + [C b3 J {G,} , (14b) 

{tom)} = [C hl l {G,} + [C b2 l {G 2 } + [C w l {G,} . (14c) 

The interface compatibility conditions will be used to evaluate {GJ, {G 2 }, and {G 3 }. Since the 
solutions of the differential equations are going to be evaluated each At, the compatibility conditions will 
be approximately satisfied by requiring the interface displacement, velocity, and acceleration of both 
bodies be equal at the end of each At. That is 


{XL(ti+ ,)} = {X [ b (t i + 1 )} ; {X' a (ti .,,)}= {X' b (ti + 1 )} ; {X'a(t m )} = {X , b Ct i + ,)} ■ (15) 

Substituting equations (5), applied to the interface coordinates, and equations (8) in equations (15) yields 

14>U {q a (tH i)} - (4>bl {Uh + .)} = r<t>b] {qfbdi ,■)} - [<!>'! Iqad, h)} , 06a) 

[cj)'a] {i(t i+ ,)} - [4>'bl {q b (t i+ .)} = M>' b ] {to, ,)} - M {toi + ,)} > 06b) 

I4»!.l {q a (ti m)} - M {ton i)} = [4> b l {q’bd, , i)} - M {to. . .)} ■ 06c) 

The terms on the left hand side of equations (16) represent the difference in displacements, velocities, 


and accelerations of the interface DOF’s of substructures a and b, due to the external applied forces. 
These terms are written as 

{8(t i+ ,)} = [<{>[,] {q a (tj + 1 )} - f<|>b] {q b (tj + i)} , (17a) 

At i+1 )} = [<J)L] {q a (tj+ i)} - [4>b] {q b (ti+ i)} , 07b) 

{5(tj , ,)} = 14>'J {q u (ti + ,)} — I4>bj {q b (t, + i)} • (17c) 
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Substituting equations (13), (14), and (17) into equations (16) gives 


{8(t l + l )} = ll4>bJlc b il - l4»aJ IQul] (Gil + lL4>hJ ic b2 ] - I4>LI IQ. 2 IJ {g 2 > 

+ [[4>L] [C b3 ] - [Vj [C a3 ]] {G 3 } , 

{5(t i+ ,)} = [M>L][c b ,i - i<j>Li [Caiii {g,} + rr«t>b] - [<i>'i [Cdi {g 2 } 

+ [[4>b] [C b3 ] - I4*a] [C a3 ]l {G 3 } , 

{5(t i+ ,)} = [[<t>' b ][C bl ] - M [C al ]] {G,} + [I4* l b ] [Cb2] - I4>'al [C a2 ]l {g 2 } 

+ [f4> b ] [C b3 ] - [<t>Ll [C a3 ]] {G 3 } . 

Equations (18) can be put in a single matrix equation in partitioned form. This results in the following 


*{8(t, , ,)} 


[C,] [C 2 ] [C 3 ] 

W 

tf(tir,)} 

= 

[C|] [&] [C 3 | 

{G 2 } 

{8(h+ 1 )} 


[C,] [C 2 ] [C 3 ] 

{g 3 }_ 


| 

- or 


- 

‘{8(t i + 

.)}“ 


’{G,f 


{S(t i + 

,)} 

= [Cl 

{G 2 } 

3 

{5(tj + 

1 )} 


{g 3 } 


{G ,}, {G 2 }, and {G 3 } can be obtained from equation (20) by inverting the coefficient matrix [C). 


(18a) 


(18b) 


(18c) 


{G,f 


{8(t i ( ,)}" 

{g 2 } 

- ter' 

{8(t 1+ ,)} 

{g 3 } 


{8(t; m)I 
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It should be observed that the coefficient matrix in equation (21) does not depend on anythin* 
related to time except the time step At. Therefore, if the At is held constant during the integration, the 

coefficient matrix used to compute {G,}, {G 2 }, and {G 3 } only needs to be computed once at the start of the 
integration. 

The essential equations for implementing the proposed transient response method have been 
developed. The task at hand now is to use these equations to form an efficient and practical algorithm for 
the solution of transient response equations. 


COMPUTATION PROCEDURE 


The computation procedure for this transient response method is best explained by providing a list 
of steps in the order they are to be implemented. The procedure consists of five initial computation steps 
that are executed once followed by eight steps that are executed once each integration time step The 
steps are as follows: 

Step I - Compute free modes and frequencies of each substructure. 

Step 2 - Select an integration time step At that is consistent with the highest substructure free 
frequency. 

Step 3 - Compute the interface compatibility matrix [C], as defined in equation (20), and its 
inverse. 

Step 4 - Set tj = integration start time, with i = I. 

Step 5 - Compute initial conditions at integration start time. These include substructure modal 
displacements and velocities, along with initial interface forces. 

Step 6 - Set tj , i = tj + At. 

Step 7 - Compute response of substructures due to applied loads at time = t ( , , by solving 
equations (9a) and (9c) with {Go} defined by equation (11). 

Step 8 - Compute interface displacements, velocities, and accelerations due to applied loads at 
time = tj f ! using equations (17). 

Step 9 - Compute the coefficients of the power series as expressed in equations (6) at time = t, , 
using equations (21). 

Step 10 - Compute the response of the substructures due to the interface forces at time = tj , , 
using equations (13) and (14). 

Step 1 1 - Compute the total response of the substructures at time = t j4 | by adding the response 
due to applied forces and the response due to interface forces using equations (8). 
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Step 12 - Compute the interface forces using equations (6) and reset initial conditions for next 
time step using equations (10) and (11). 

Step 13 - Set t, = t i+ , and return to step 6. 

Care must be exercised in carrying out step 1 to prevent an unacceptable loss of accuracy . If on< - 
simply computes the free modes and frequencies from the discrete coordinate substructure mode s an 
truncates the number of modes based on a cut-off frequency, the motion and forces at the interface wil 
not be represented very well. This is the same problem that occurs when free modes are used in modal 
synthesis methods. This loss of accuracy can be eliminated by introducing an additional step. The step 
consists of transforming the discrete coordinate models to the Craig-Bampton [6] form using the interface 
DOF’s as boundary DOF’s. The modal truncation is introduced in the formation of the Craig-Bampton 
model by truncating the fixed boundary modes based on a cut-off frequency that is consistent with the 
applied forcing function. The free modes and frequencies are computed from the Craig-Bampton models 
and are not truncated. This procedure assures the motion and forces at the interface are accurate) 

represented. 


DEMONSTRATION PROBLEM 


To demonstrate the transient response method for linear coupled substructures, a cantilevered 

2XS B) of'eTuat ZaZnO !" « 

szazz 

strained normal modes were retained up to 100 Hz. Interlace uur s were iciai 

free modes and frequencies of the substructures were computed by an eigenvalue ana y^ . * 

R minton models This eigenvalue analysis produced frequencies in the range of 4 t - > 

anTS io * Hz for My B. All modes from foe eigenvalue analysis were returned for this study. 

A benchmark model composed of the two reduced substructures was developed to investigate the 
coupled substructure transient response method. The benchmark model was formulated by !"S the 
two^ reduced beam substructures together using the direct stiffness method. The coupled equation 
motion for the benchmark model were then uncoupled by an eigenvalue analysis. The uncoup e 
equations of motion resulted in a frequency range of 1 to 185 Hz, with all modes retaine . 

Damping is normally introduced in the uncoupled equations of motion in terms of a percent of 
critical damping For this study, a value of zero percent of critical damping was used in order to compare 
the benchmark closed- form solution to the coupled substructure solution. 

To solve the transient response analysis of the two-body differential equations, a Fortran corn 
puter routine was written utilizing the computational procedure described in the previous section. For t e 
benchmark differential equations, a standard closed-form transient response routine was used. The 
applied force defined in Figure 2b was linear interpolated for both the closed-form ana ysis an 
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E I - 52080 lbs-in 2 2 2 

m - .00307169 Ibs-sec. /in. 

L - 48.0 in. 

Figure 2a. Cantilever beam with applied load F(t). 


F(t) 



Figure 2b. Definition of applied force. 


Body A Body B 

Figure 2c. Two-dimensional finite element model of cantilevered beam. 
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method presented. The highest frequency of the two substructures is 286 Hz. This gives a minimum 
period of 0 0035 seconds. In order to achieve accurate results in transient analysis, a time step must be 
chosen below the minimum period of the problem being solved. For this study, a range of time steps was 
chosen to research the accuracy of the proposed method. The time steps chosen range from 0.0003 to 
0 0042 seconds. The benchmark transient response was computed using the same time steps used in the 
substructures transient response method. All transient response analyses were run from 0 to 1 .0 seconds. 

Interface shear force and moment were recovered from the analysis of the two substructures and 
were compared to the benchmark results. Also compared were the tip acceleration, velocity, and dis- 
placement Using a time step of 0.001 seconds, the tip acceleration, velocity, and displacement, obtained 
from the response analysis of the coupled substructures and the benchmark response analysis were 
plotted versus time in Figures 3a, 3b, and 3c. Within the fidelity of the plots, the results fell on top of 

each other. 

The solutions from the benchmark closed-form transient response are compared to the substruc- 
ture method by the relative error defined in the following equation: 


RMS Error = 


Big 


2 (X b -X s ) 2 

NT 


NT 


where X b is the benchmark time vector result (i.e. , force, acceleration, etc.); X s is the substructures time 
vector result (i.e., force, acceleration, etc.); NT is the total number of time points computed; and Big is 
the absolute largest benchmark result. 

The RMS error results are plotted versus the time step used in the transient response analysis. The 
RMS errors for the tip acceleration, velocity, and displacement results are shown in Ftgure 4a, RMS 
errors for the interface shear force and moment are shown in Figure 4b. A corre at, on extsts between the 
order of the lime derivative results to that of the RMS error in Figures 4a and 4b. For example the Up 
acceleration is an order of time derivative higher than the tip velocity, and the RMS error of e ip 
acceleration is an order of magnitude higher than that of the tip velocity. Ltkewtse, the tip velocity RMS 
error is higher compared to the tip displacement. This can also be seen between the interface shear force 
and moment where an order of one exists between the time derivatives. 

For a reasonable time step, the coupled transient response method presented gives accurate results 
of the demonstration program. As the time step gets closer to the minimum period, a leveling off of he 
RMS error occurs. This effect corresponds to accuracy of the integration for a given time step used 
analysis. Using a time step larger than the minimum period results in erroneous results. 
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Figure 3a. Ti 






1.00E+00 



Figure 4a. RMS error versus time step. 
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Figure 4b. RMS error versus time step. 


SUMMARY 


A method for transient response analysis of linear coupled substructures has been presented. The 
equations of motion are solved in a stepwise fashion with the interface forces between substructures 
being approximated by a third-order power series in time. The coefficients in the power series are 
evaluated for each time step by enforcing compatibility of displacements, velocities, and accelerations at 
the substructure interfaces at the end of each time step. The method is approximate because the interface 
compatibility is enforced only at the end points of a time step rather than continuously. Therefore, as in 
most numerical integration methods, the smaller the time step the more accurate the results. The evalua- 
tion of the coefficients in the interface force power series is accomplished by a matrix multiplication. The 
only dependence on time in the interface compatibility coefficient matrix is the integration time step. 
Therefore, the matrix is generated only once unless the time step size is changed. The accuracy of the 
method has been validated by comparing the results with a closed-form response analysis of cantilever 
beam. 


The interface compatibility between substructures is maintained by the use of the interface 
compatibility coefficient matrix. Therefore, this method is especially appealing for the response analysis 
of structural systems that have substructure interface boundary conditions that change during the 
response run. The coefficient matrix is simply changed or regenerated to reflect the change in boundary 
conditions between the substructures. The application of this method to the response analysis of a launch 
vehicle while lifting off from the launch pad was mentioned earlier. Another application is for structures 
that exhibit slip-stick motion at interfaces caused by friction. The appropriate interface compatibility 
matrix is selected by monitoring the interface force and determining if the interface is slipping or 
sticking. 

Since the equations of motion of each substructure are solved independently for each time step, it 
appears the method would be ideally suited for parallel processing. 
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